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I. INTRODUCTION 

Graphene, a single sheet of carbon atoms in a hon- 
eycomb lattice has attracted enormous interest since 
its experimental discovery and isolation. 1,2 Its quasi- 
particle states behave like massless "relativistic" Dirac 
fermions^ and obey a two-dimensional Dirac equation, 
where the light speed is replaced by the Fermi veloc- 
ity that is 300 times smaller than speed of light in vac- 
uum. Beside various applications in electronic devices, 
graphene physics opens wide research field unifying low- 
energy condensed matter physics and quantum electro- 
dynamics (QED)^r— Many fundamental nonlinear QED 
processes, specifically, electron-positron pair production 
in superstrong laser fields of ultrarelativistic intensities, 10 
observation of which is problematic yet even in the cur- 
rent superintense laser fields, have their counterparts 
in graphene where considerably weaker electromagnetic 
fields are required for experimental realization of "anti- 
matter" production from vacuum. In this connection one 
can note Klein paradox^ - — Schwinger mechanism-^— 
and Zitterbewegungi^r— for particle-hole excitation, as 
well as diverse physical and applied effects based on Zit- 
terbewegung, e.g., minimal conductivity at vanishing car- 
rier concentration 2 ^^ etc. 

Due to massless energy spectrum the Compton wave- 
length for graphene quasiparticle tends to infinity. On 
the other hand, in QED the Compton wavelength is char- 
acteristic length for particle-antiparticle pair creation 
and annihilation. So, at the interaction of an electro- 
magnetic field with an intrinsic graphene there is not a 
quasiclassical limit, since no matter how weak the applied 
field is and how small the photon energy is, the particle- 
hole pairs will be created during the whole interaction 
process -at arbitrary distances. 

In graphene wave-particle interaction can be charac- 
terized by the dimensionless parameter 

eEvp 1 

which represents the work of the wave electric field E 
on a period 1/uo in the units of photon energy e 7 = Hoj. 
Here vf is the Fermi velocity: ~ c/300, and e is the 



elementary charge. The average intensity of the wave 
expressed by x, can be estimated as: 

I x = x 2 x 3.07 x 10 11 W cn:T 2 [/iw/eV] 4 . 

Depending on the value of this parameter one can 
distinguish three different regimes in the wave-particle 
interaction process. Thus, x << 1 corresponds to one- 
photon interaction regime, x ~ 1 - to multiphoton in- 
teraction regime, and x >> 1 corresponds to static 
field limit or Schwinger regime. As is seen, the inten- 
sity I x strongly depends on the photon energy. Par- 
ticularly, for infrared photons: e 7 ~ 0.1 eV, multipho- 
ton interaction regime can be achieved at the intensities 
I x — 3.07 x 10 7 W cm~ 2 . Note that in case of free elec- 
trons at the same photon energies multiphoton effects 
take place at the intensities I ~ 10 16 W cm~ 2 i - Such a 
huge difference, as well as the gapless particle-hole energy 
spectrum in graphene makes realistic another interesting 
nonlinear QED process: multiphoton excitation of Dirac 
vacuum with Rabi oscillations at ultrafast excitation^ 
In the present work the creation of particle-hole cou- 
pled states in graphene via multiphoton resonant excita- 
tion by laser fields is studied. It is well known that Rabi 
oscillation of states' populations is the coherent response 
of two-level atomic systems under resonant excitation. 
The one-photon resonant excitation of atoms and associ- 
ated Rabi oscillations have been comprehensively studied 
both theoretically and experimentally and described in 
numerous review articles and books (see, e.g.j^l). Similar 
phenomena have also been observed in semiconductors. 25 
Recently, the Rabi oscillations in graphene at one-photon 
interband excitation (at x << 1) and its influence on 
the dynamic conductivity were investigated in Refs. [26l — 
[28| . On the other hand, at a wave-particle interaction in 
graphene due to free-free intraband transitions we have 
a situation analogous to resonant excitation of quantum 
systems with permanent dipole moments, where direct 
multiphoton transitions are very effective i22r— Hence, it 
is of interest the creation of particle-hole coupled states 
in graphene via multiphoton resonant excitation. This 
process apart from fundamental interest may also have 
practical applications. In particular, particle-hole an- 
nihilation from coherent superposition states will cause 



2 



intense coherent radiation of harmonics of the applied 
wave-field. We consider multiphoton interaction regime 
at x ~ 1- Accordingly, the time evolution of consid- 
ered process is found using a nonperturbative resonant 
approach arising from the quantum kinetic equations. 

The paper is organized as follows. In Sec. II the set of 
equations for single-particle density matrix is formulated. 
In Sec. Ill we present the solution of those equations in 
the multiphoton resonant approximation. In Sec. IV 
the results of corresponding numerical simulations are 
presented. Finally conclusions are given in Sec. V. 



We write the fermionic field operator in the form of an 
expansion in the free Dirac states: 

*(x,y,t) =^a p , CT (i)*p >CT (a;,y), (4) 

p,<T 

where the creation and annihilation operators, a (t) 
anda Pi(T (i), associated with positive (a = 1) and negative 
(a = —1) energy solutions satisfy the anticommutation 
rules at equal times 

{a P)CT (*)>apv(0}i=t' =<W"W , (5) 



II. BASIC MODEL AND EVOLUTIONARY 
EQUATION FOR SINGLE-PARTICLE DENSITY 
MATRIX 



Let graphene interacts with a plane quasimonochro- 
matic laser radiation of carrier frequency ui and slowly 
varying envelope. To clarify the picture of stated in 
this paper problem for creation of coherent superposition 
states, we consider the case of interaction when the laser 
wave propagates in perpendicular direction to graphene 
plane (XY) to exclude the effect of magnetic field. This 
travelling wave for electrons in graphene becomes a ho- 
mogeneous time-periodic electric field. Let it to be di- 
rected along the X-axis with the form (constant phase 
connected with the position of the wave pulse maximum 
with respect to graphene plane is set zero): 



E (t) = xE cos uit. 



(1) 



The problem that we attempt to solve in the given field 
approximation is analogous to Rabi oscillations in two 
level atomic systems with permanent dipole moments. In 
this case, the physical picture of resonant wave-graphene 
interaction will be more transparent in the length gauge. 
So, for the interaction Hamiltonian we will use length 
gauge describing the interaction by the potential energy. 
Cast in the second quantization formalism, the Hamilto- 
nian is 



H = / ^ + H s ^dxdy, 



(2) 



where "J/ is the fermionic field operator, H s is the single- 
particle Hamiltonian in the external homogeneous elec- 
tric field ([T]) . Omitting here real spin and valley quantum 
numbers, single-particle Hamiltonian can be written as: 



v F 



Px~ Wy 

Px + Wy 



exE 
exE 



, (3) 



where vp ps c/300 is the Fermi velocity (c is the light 
speed in vacuum), p = {p x ,p y } is the electron momen- 
tum operator. The first term in Q is the Hamiltonian of 
two-dimensional massless Dirac fermion and the second 
term is the interaction part. 



{2p )CT (*),a p , i(T , = {a P At),a p/ ^(t')} t=t , =0. 

The free Dirac solutions corresponding to energies £ a 



(6) 

: solutions corresponding to energies £, 
a ^F\jp 2 x +P 2 y (<7 = ±1) are: 

* P A*,y) = ^( ae m P) )e*^ + ™\ (7) 

where S is the quantization area (graphene layer surface 
area) and 



6 (p) = arctan ( — 

\Px 



(8) 



is the angle in momentum space. 

Taking into account Eqs. (HJ-®, the second quantized 
Hamiltonian can be expressed in the form: 



H = E ° (P) «p<t«p^ 



^W^E D °°' (P'P')«p )CT V,a'. (9) 



where 



Daa> (P,p') = — 



1 + aa'e 



i(e( P ')-e(p)) 



(10) 



We will use Heisenberg representation, where operators 
evolution are given by the following equation: 



ih 



dL 
~dt 



L, H 



(11) 



and expectation values are determined by the initial den- 
sity matrix D: 



<L>= 



= S P ( 



DL 



(12) 



The single-particle density matrix in momentum space is 
defined as: 



/W 2 (pi,P2,t) =< a+ >0 . 3 (t)op llffl (f) > 



(13) 
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For initial state of graphene quasiparticles we assume 
ideal Fermi gas in equilibrium. This means that the ini- 
tial single-particle density matrix is diagonal and we have 
Fermi-Dirac distribution: 



/w(p,p',0) 



1 + e t 



■-. ■■■■ -Jj 'Vl>'*""' 



(14) 



Here p is the chemical potential, T is the temperature in 
energy units. Taking into account definition (fT3"]) , from 
Eq. (|TT|) one can obtain evolution equation for single- 
particle density matrix: 

■ t -dp cricr2 (pi,p2,t) 

™ = (Pl) ~ £a 2 (P2)\ Po t a 2 (Pl,P2, t) 



dt 



p. a 



-D aia (Pl,p) Paa 2 (P,P2,t)} . 



(15) 



Then taking into account the following relation with 
Dirac delta function S (a) : 



and substitution 



xe 



E 



d 

' dx = 2iri—5 (a) 
da 



S 



dp, 



we obtain closed set of equation for the density matrix 
elements: 



-eE(t) 



dt 

dp<T,a(P,P2,t) 



Pl=P 



2x 



eE (t) dO (p) 



dp x 



x W-<t(p, P, t) - P- CT)CT (p, P, t)] , 

9p CT ,-<r(p, P, t) _ ^ ,^ 9/9 CT) _ ff (pi,p,t) 



(16) 



-e£(i) 



<9t 

5p<r-<r(P,P2,*) 



dpix 



pi=p 



dp2x 

eE (t) dO (p 

2i dp x 



= -r£r {p)pa-a(P,P,t) 

in 



[p<tAp> P> *) - -o-(p, P, *)] • ( 17 ) 



In Eqs. HH]) and (JT7J) one can eliminate the terms with 
partial derivatives d/dp x by the method of characteris- 
tics. The characteristic of these equations is the classical 
equation of motion: 



dp 
~d~t 



-eE (t) 



with the solution 

Px = POx +PE(t); Py = P()y, (18) 

where 

p E (t) =-e E (t') dt' = -sinivt 

Jo ^ 

is the momentum given by the wave-field. Thus, with the 
new variables Pq x , po y , and t Eqs. (1161) and (I17|) read: 



dpa,a{PQ,Po,t) = £ 

dt ~ 2 



F{po,t) 



[p<j-*(Po, Po, t) - p_ CT , CT (p , po, t)} , (19) 



5p CT ,- CT (po,Po,t) 2 ~ 

^ = -ffa (PO, f) Pa ,-<t(P0, PO, t) 



+ 7; F (Po> t) [p<r,a(Po, Po, *) - P-a -<x(po, Po, *)] , (20) 



where 



and 



F(po,t) = - , (21) 

(Pox+PeW) +Po„ 



Mpo,*) =av F ^(p x+PE(t)) 2 +p 2 0y . (22) 
Taking into account Eq. (p~5]) . it is easy to see that 



P CT ,<r'(Po,Po,0) = p CTiO ./(p,p,0). 



(23) 



To be more precise in the set of equations ([T6"l) and (fl7|) 
one should add the terms describing relaxation processes. 
Since we have not taken into account the relaxation pro- 
cesses, this consideration is correct only for the times 
t < T m i n , where r m i n is the minimum of all relaxation 
times. Thus, full dynamics in the absence of any losses 
is now governed by Eqs. (|19[) and ([20)) . These equations 
yield the conservation law for the particle number: 

Pi,i(Po, Po, t) + P-i -i(Po, Po, t) = 



= Pi.i(Po, Po, 0)+p_i -i(p , Po, 0) = 5 (p , p, T) . (24) 

Here we have introduced the notation 5 POjAIj t, which ac- 
cording to Eq. ffl4|) is: 



1 + e t 



1 + e" 



In Eqs. (fT^j) and (|2U)) diagonal elements represent 
particle Af (po,t) = pi,i(po, Po> t) an d h°l e (1 — 
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/ o_i l _i(po,po,i)) distribution functions, while nondiag- 
onal elements Pi,-i(po, Poj t) = P-% i(Po>Po>*) describe 
particle-hole coherent transitions. Introducing the new 
quantity (J (po,*)): 



Pi,-i(PO)Po,*) = ij(p ,t)exp\ -i- / £i(p ,t')dt' 



and taking into account that p_i ) _i(po, po, t) = S Po . m .t — 
Af (po, i), from Eqs. ([T9|) and (|20|) we obtain the following 
set of equations: 

gAOpo^ _ 1 

a* --2 F(P0 '* } 



XitlbJ 



J (p 0> t)exp{-i- J £ x (p ,t')dt' } + < .< . 



&7(Po,t) 1 



. (25) 



dt 



F{p ,t)exp\i- / fi(p ,f)^ 



x (p 0) f) - Ej, ,u,t] 



(26) 



This set of equations should be solved with the initial 
conditions: 



J(po,0)=0; AA(p o ,0) 



1 + e T 



(27) 



III. MULTIPHOTON RESONANT EXCITATION 

The Eqs. (|2"5")) and (j2l)]) represent linear set of equations 
with periodic coefficients and are analogous to Bloch 
equations^ which describe Rabi oscillation of states' 
populations of two-level atomic system under resonant 
excitation. Note that there are significant differences be- 
tween usual Bloch equations and Eqs. (|2"5|) . (|2l)|) . Thus, 
from the Floquet theorem and Eq. (|22|) follows that in- 
stead of stationary levels vppg and —vfPo due to free- free 
intraband transitions we have quasistationary states with 
quasienergies: £± (s) — ±£e (Po) + sw ; s = 0,±1,±2, .... 
Here 



£e (Po) = ^ 



2tt/u 



£\ (Po.*) dt 



Vp 

2tt 



P ° X ~ ~~ sinUjt ) +P0ydt l-' S ) 




FIG. 1: Graphene conical dispersion with interband multi- 
photon transitions induced by external electric field. 



F (po, t + 2tt/oj) = F(p ,t), which in tern contains all 
harmonics of driving field. In usual Bloch equations cou- 
pling contains only fundamental oscillations, which pro- 
vide only direct one photon resonant excitation. Situa- 
tion is analogous to resonant excitation of systems with 
permanent dipole moments, where, as has been shown 
in Ref. [29(, it is possible to decouple slow and rapid 
oscillations and to disclose the resonant dynamics of the 
wave-particle interaction. 

Because of space homogeneity of the field (pj, the 
generalized momentum of a particle conserves, so that 
the real transitions in the field occur from a —£e (Po) 
negative energy level to the positive +£e (po) energy 
level and, consequently, the multiphoton probabilities of 
particle-hole pair production will have maximal values 
for the resonant transitions 



2£ Eo 



xhjj. 



(29) 



Note that the resonant condition (|29|) is equivalent to 
crossing of Floquet ladders: £_ (s + n) ~ £ + (s). Figure 
1 schematically illustrates multiphoton interband transi- 
tion between the two states in the filled lower cone and 
the empty part of upper cone. To decouple slow and 
rapid oscillations at the resonant condition (|2"9")l . we fol- 
low the ansatz developed in Ref. [1^] and represent Eqs. 
(1251) and M in the form: 



dAf(pg,t) 
dt 



x G *s (Po. Ea) e-Kt^Cpo)-^)* + c c ; (30 ) 



is the mean value of classical energy in the field ([T}, 
which has a nonlinear dependence on the amplitude of 
the wave-field. The latter in the physical sense is the 
dynamic Stark shift due to free-free intraband transi- 
tions. Two Floquet ladders are coupled by the term 



dt 2 

S 

x [2Af (p ,f) - S P0!jU , T ] 



(31) 
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Here s photon coupling coefficient is: 



G s (p ,E ) = 5- / F(p , 
2tt Jo 



t) 



x exp 



m#0 



2£m (Po) -imut ( ,.^-1 
ihlO 



dt. 



(32) 



where 



£ m (Po) 



2tt 



2tt/u 



£ 1 (p ,t)e imut dt. (33) 



Close to resonance (|2"9"|) . the main coupling term in Eqs. 
(130)) and (|31[) becomes slowly varying term with s = n. 
The remaining nonresonant and rapidly oscillating terms 
are responsible only for dynamic Stark shifts^ Thus, for 
time average functions Af (po,t) and J (po,t) one can 
obtain the following set of equations: 

m ^ t) = -ij(p ,t)G; i (p ,i?o) e ^" t + c.c., (34) 



where 



A 



1 



1 + e t 1 + e 

is the initial population inversion, and 



(40) 



yJ\G n (po,E )\ 2 + (6 n + 6 st ) 2 (41) 



is the generalized Rabi frequency. The solution (j3"9")l ex- 
presses Rabi flopping among the particle-hole states at 
the multiphoton resonance. The solutions (|38|) . (|39| have 
been derived using the assumption that A/"(po,i) and 
J' (po, t) are slowly varying functions on the scale of the 
wave period, which put the following restrictions: 



(|G n ( Po ,£7o)|, W, \Sst\) «u 



(42) 



on the characteristic parameters of the system consid- 
ered. 

For the exact resonance (5 n + 6 s t = 0) we have f2 n = 
\G n (po,Eq)\ and the solutions become: 



dj (p ,t) 
fit 



+ i<5 s4l 7(po,t) 



= l -G n (p ,So)e i5 "' (2AT( Po ,i) - S Po , M , T ) 

Here we have introduced resonance detuning 

2£e (po) - nhu 
o n = 7 



and dynamic Stark shift 
5 st 



1 ^|G s (p ,^o) 



2uj (n — s) 



(35) 



(36) 



(37) 



The latter is the result of virtual nonresonant transition 
within Floquet states. Thus, we have a set of linear or- 
dinary differential equations the solution of which at the 
initial condition (l27l) is: 



-7f, n .\ iS n t ^n (P0>ffi>) A 

J (p ,t) = e — A POiM ,t 



■ f sin <),,/ - i — ^ ^ st (1 — cosf2„t) I , (38) 



TTf n f\- S P0-M,T , \Gn (PO,^o)| a 

Al (Po,i) - — ^ 1 ^ A P0,M,T 



\G n (P0,^0) 



cos f2„ i 



(39) 



Ar(po,t) = ^^cosO n i+^^. (44) 

At one-photon interband excitation (when %o << 1) (|3"2")l 
one can omit nonlinear over _Eo terms and for Rabi fre- 
quency we have: 



9,1 



eE | sin 9 (p )| 
2p 



Taking into account resonant condition 2poVp ~ fiw, the 
latter can be expressed through the interaction parame- 
ter xo = eE v F /(haj 2 ): 

fli =w|sine(po)|xo- 

With increasing of pump wave intensity, the Rabi oscil- 
lations appear corresponding to multiphoton transitions. 
At that one should take into account the intensity effect 
of the pump wave on the quasienergy spectrum (Stark 
shift due to free-free intraband transitions) and the dy- 
namic Stark shift due to virtual nonresonant transitions. 
For x ~ lj the probabilities of multiphoton transitions 
are essential up to photon numbers n ~ 5. For this pho- 
ton numbers the Stark shift (|38|) is not essential, while 
the modification in the quasienergy spectrum is consid- 
erable. Isolines of quasienergy spectrum, defined by Eq. 
(|2"5|). are no longer circles but ellipse-like. In Fig. 2, iso- 
lines of quasienergy corresponding to resonant condition 
2£ Eo (p ) /hw « n flU with detuning \6 n \/u = 0.02, 
for n — 2,3,4 are shown in the multiphoton interaction 
regime: xo — 1- As is seen, modification of unperturbed 
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FIG. 2: Isolines of quasienergy corresponding to resonant con- 
dition 2£ Eq (po) /tiu ~ n (J29J) with detuning |<5„| /u = 0.02, 
for n = 2, 3, 4. The electric field dimensionless parameter is 
taken to be xo = 1.0. The momentum components are nor- 
malized to hw/vF- 



energy spectrum in the direction perpendicular to the 
electric field vector, is considerable (for example, at n = 3 
we have 13% deviation). Thus, in the multiphoton inter- 
action regime one should expect photoexcitation of parti- 
cle distribution function just along the modified isolines, 
in accordance with Eq. ((44)) . 



IV. NUMERICAL TREATMENT 

We have also performed numerical simulations and in- 
tegrated Eqs. (piS]) and with the fourth-order adap- 
tive Runge-Kutta method. As we mainly interested by 
interband multiphoton transitions, for all calculations the 
chemical potential and temperature are fixed and are 
taken to be fi/hu = 0.1 and T/fku — 5 x 10~ 3 . 

In Figures 3-5 photoexcitations of Fermi-Dirac sea is 
presented: density plot of the particle distribution func- 
tion Af (po, t) is shown for various pump wave intensities 
and instants. In Fig. 3 corresponding to \o — 0.02, we 
see only creation of particle-hole pair in graphene at the 
one-photon resonant excitation. In Fig. 4, the pump 
wave intensity is larger: xo — 0.5 and, as a consequence, 
we see resonant rings corresponding to multiphoton exci- 
tation up to four photons. In Fig. 5, which corresponds 
to Xo = 1) it is clearly seen also the ring for five photon 
resonance. At that, the ring for one-photon excitation is 
smeared, since Stark shift for this energy is comparable to 
hcu and the condition (|42l) for resonant Rabi oscillations 
at one-photon excitation is not fulfilled. As is seen from 
Figs. 4 and 5, the excitation of Fermi-Dirac sea takes 
place along the ellipse-like isolines of quasienergy spec- 
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FIG. 3: (Color online) Creation of particle-hole pair in 
graphene at the one-photon resonant excitation. Particle dis- 
tribution function M (po , t) (in arbitrary units) at instant 
U = 25T, as a function of scaled dimensionless momentum 
components {poxVF /Hu),poyVF /Hu>}. The electric field dimen- 
sionless parameter is \o = 0.02. 
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FIG. 4: (Color online) Creation of particle-hole pair in 
graphene at the multiphoton resonant excitation. Particle 
distribution function JV"(po,£) (in arbitrary units) at instant 
U = 100T, as a function of scaled dimensionless momentum 
components {po x v f / hw , po y v f / fi^} ■ The electric field dimen- 
sionless parameter is xo = 0.5. 
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FIG. 7: (Color online) Same as Fig. 6, but for three-photon 
resonance (n = 3). The angle is taken to be O (po) = 7r/2 rad. 
The electric field dimensionless parameter is \o = 1-0. 



FIG. 5: (Color online) Same as Fig. 4, but for stronger wave- 
field with \o = !• 
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FIG. 6: (Color online) Two-photon resonance (n = 2). Rabi 
oscillations of the particle distribution function M (po,t) for 
the fixed angle 6 (po) =0.2 rad versus the scaled dimen- 
sionless momentum povf/&^- The electric field dimensionless 
parameter is \o = 0.5. 



trum defined by Eq. (I28I) . in accordance with analytical 
treatment (see, Fig. 2). 

To show the dynamics of multiphoton excitation of 
Fermi-Dirac sea in Figs. 6-8, we present Rabi oscilla- 
tions of the particle distribution function N (po, t) for the 
fixed angles versus the scaled dimensionless momentum 
Povp/huj. Figure 6 corresponds to two-photon resonance 
for the angle (p ) =0.2 rad and xo = 0.5. It is clearly 
seen Rabi oscillations of J\f (po,t) with the mean period 
Tr = 22T (T is the wave period). Rabi oscillations for 
three-photon resonance is displayed in Fig. 7 for the 
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FIG. 8: (Color online) Four-photon resonance (n — 4). Same 
as Fig. 7, but for the angle G (po) = 0.6 rad. The electric 
field dimensionless parameter is xo = 1.0. 



angle (po) = 7r/2 rad and xo = 1- Here mean Rabi 
period is Tr = 157". Four- photon resonant Rabi oscil- 
lations with the mean period Tr = 28T, at the angle 
(po) = 0.6 rad is shown in Fig. 8. 

To show the dependence of Rabi oscillations on the 
angle 0(po), in Figs. 9 and 10 colored 4D density 
plot of Rabi oscillations of the particle resonant distribu- 
tion function Af r (po, t) on isosurfaces 2£e (Po) = n, 
for n — 3,4 are shown at xo = !• As is seen from 
these figures, for odd photon resonance the excited dis- 
tribution function Af r (po , t) is maximal at the angle 
(po) — tt/2 rad (perpendicular to applied electric field 
vector), while for even photon resonance at (po) = 
7r/2 rad we have a minimum and main excitation takes 
place close to (po) = tt/4. The latter is connected 
with the fact that coupling term in Eqs. (|3TJ|) . (|3"Tj) at 
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JVr(Po.f) 




FIG. 9: (Color online) Colored 4D plot of Rabi oscillations of 
the particle distribution function Af r (po,t) for three-photon 
resonance on isosurface 2£e (po) /fow = 3. The electric field 
dimensionless parameter is xo = 1-0. 




FIG. 10: (Color online) Same as Fig. 9, but for four-photon 
resonance: 2£e (po) /fku = 4. 
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(po) = 7r/2 rad contains only odd harmonics of the 
pump wave: G 2s (0,p Oyi E o ) = 0. 

Summarizing, we see that numerical simulations are 
in agreement with analytical treatment in multiphoton 
resonant approximation and confirm simple physical pic- 
ture described in previous section. We see that with 
laser fields of moderate intensities one can observe mul- 
tiphoton excitation of Fermi-Dirac sea, the probabilities 
of which are comparable to and sometimes larger than 
one-photon excitation probability. 

V. CONCLUSION 

We have presented a theoretical treatment of the co- 
herent nonlinear response of a graphene under multipho- 
ton interband excitation by a laser radiation. The evo- 
lutionary equation for single-particle density matrix is 
formulated arising from the second quantized formalism. 
The time-dependent single-particle density matrix in the 
given field of a laser radiation is calculated in the multi- 
photon resonant approximation. The Rabi oscillations of 
Fermi-Dirac sea at multiphoton excitation depending on 
the time, momentum, and photon number have been con- 
sidered and analyzed also on the base of numerical simu- 
lations. The results obtained demonstrate that the Rabi 
oscillations of Fermi-Dirac sea corresponding to multi- 
photon excitation can already be observed for such laser 
fields where the work of electric field on the wave period is 
comparable to photon energy e 7 = fvjj. For mid- infrared 
lasers e 1 ~ 0.1 eV, multiphoton interaction regime can 
be achieved at the intensities I > 10 7 W cm~ 2 , for the 
time scales 1.0 ps. For near-infrared range of frequen- 
cies Ey ~ 1. eV, multiphoton interaction regime can be 
achieved at the intensities / > 10 11 W cm~ 2 , for the time 
scales 100.0 fs. 
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